BNIPL is a promising biomarker of laryngeal cancer: novel insights from bioinformatics analysis and experimental validation

Background Laryngeal cancer (LC) is a malignant tumor with high incidence and mortality. We aim to explore key genes as novel biomarkers to find potential target of LC in clinic diagnosis and treatment. Methods We retrieved GSE143224 and GSE84957 datasets from the Gene Expression Omnibus database to screen the differentially expressed genes (DEGs). Hub genes were identified from protein-protein interaction networks and further determined using receiver operating characteristic curves and principal component analysis. The expression of hub gene was verified by quantitative real time polymerase chain reaction. The transfection efficiency of BCL2 interacting protein like (BNIPL) was measured by western blot. Proliferation, migration, and invasion abilities were detected by Cell Counting Kit-8, wound-healing, and transwell assays, respectively. Results Total 96 overlapping DEGs were screened out from GSE143224 and GSE84957 datasets. Six hub genes (BNIPL, KRT4, IGFBP3, MMP10, MMP3, and TGFBI) were identified from PPI network. BNIPL was selected as the target gene. The receiver operating characteristic curves of BNIPL suggested that the false positive rate was 18.5% and the true positive rate was 81.5%, showing high predictive values for LC. The expression level of BNIPL was downregulated in TU212 and TU686 cells. Additionally, overexpression of BNIPL suppressed the proliferation, migration, and invasion of TU212 and TU686 cells. Conclusion BNIPL is a novel gene signature involved in LC progression, which exerts an inhibitory effect on LC development. These findings provide a novel insight into the pathogenesis of LC. Supplementary Information The online version contains supplementary material available at 10.1186/s12920-024-01811-z.


Introduction
Laryngeal cancer (LC) is a kind of malignancy of head and neck with high incidence, with approximately 1.7 million cases of LC and nearly 90,000 deaths worldwide each year [1].More than 40% of LC patients with LC are diagnosed at an advanced stage, often with spread to the supraglottic and para glottic space [2,3].The 5-year overall survival rate for LC shows continued decline in the last few decades [4].The occurrence of LC is multifactorial and is related to common carcinogens (tobacco and alcohol) [5].Radiotherapy, chemotherapy, and surgical resection are the main treatments for LC, but the recurrence rate of advanced LC is 25-50%, and there are still numerous patients with poor prognosis [6,7].Therefore, screening the promising biomarkers and therapeutic targets for the early diagnosis of LC has become a hot topic.
Currently, the potential for genomic data sharing plays a crucial role in cancer research [8].The identification of genes involved in cancer progression holds great promise for medical diagnosis [9].After retrieving quantitative and clinical gene expression data of LC from the Cancer Genome Atlas and Gene Expression Omnibus (GEO) database, the role of B7-H3 in the diagnosis and prognosis of LC is determined [10].Based on the robust rank aggregation analysis for the consolidation of gene expression datasets from the GEO database for patients with LC, CDK1, PC24, HOXB7, and SELENBP1 are screened as potential biomarkers of LC [11].Therefore, searching for potential gene markers as the therapeutic targets for LC is utmost essential.
BCL2 interacting protein like (BNIPL) is a pro-apoptotic gene that can maintain homeostasis between cell proliferation and apoptosis in vivo [12].There are 2 spliceosomes of BNIPL, BNIPL-1 and BNIPL-2.Overexpression of BNIPL-2 can significantly inhibit tumor growth by regulating molecules related to cell proliferation and apoptosis of hepatocellular carcinoma cells [13].Additionally, BNIPL-2 can interact with Bcl-2 and Cdc42GAP to regulate cell apoptosis in hepatocellular carcinoma [14].However, whether BNIPL is involved in the progression of LC remains unknown and need to be further explored.
In this study, we screened out common differentially expressed genes (DEGs) from the GSE143224 and GSE84957 datasets.The hub genes of LC were determined after constructing protein-protein interaction (PPI) network.Then we investigated the role of BNIPL in proliferation, migration, and invasion of LC in vitro.It is hoped that this study will provide novel insights of LC pathogenesis and BNIPL can be a potential target of LC in clinic diagnosis and treatment.

Microarray dataset selection
The gene expression datasets of LC-related genes were downloaded from GEO of National Center for Biotechnology Information (https://www.ncbi.nih.gov/geo/).We searched "Laryngeal Squamous Cell Carcinoma" in the GEO database, and the LC-related microarray datasets, GSE143224 (11 control samples vs. 14 tumor samples) and GSE84957 (9 control samples vs. 9 tumors samples) were selected for analysis of this study.

Identification of DEGs
GSE143224 and GSE84957 datasets were analyzed using Gene Expression Profiling Interactive Analysis (GEO2R; www.ncbi.nlm.nih.gov/geo/geo2r).|logFC|≥ 2 and adj.p ≤ 0.05 were set as the cutoff criteria for selecting DEGs.The identified DEGs were visualized using heat maps and volcano maps, and the samples were normalized and corrected by boxplot.

Enrichment analysis of DEGs
DAVID (The Database for Annotation, Visualization and Integrated Discovery) (http://david.abcc.ncifcrf.gov/)database was applied to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.GO analysis results were divided into molecular function (MF), biological processes (BP), and cellular component (CC).The results with the minimum p value were regarded as the most significantly enriched targets, which were displayed using the enrichment analysis bar chart and enrichment analysis bubble chart.

PPI network construction and hub genes identification
PPI network were constructed based on STRING online database (https://www.string-db.org/) to explore the interaction network between proteins encoded by DEGs.The confidence interaction score was set at 0.4 as the significance criterion.Visualization of the protein interaction network was carried out using Cytoscape software.The significant modules were screened out from PPI networks using Molecular Complex Detection online tool.The plugin CytoHubba (Version 0.1) of Cytoscape was used to calculate the degree of protein node, and the hub genes were selected according to the connection degree.

Analysis of hub genes
GO chord plots of hub genes were draw to reveal the relationship between proteins and the changes of the functional pathways.The expression of hub genes in GSE84957 dataset was used as a variable to perform principal component analysis (PCA).Then the expression levels of hub genes in GSE84957 dataset were visualized using ridgeline plots.Above analyis and visualization were performed using R package ggplot2.
The receiver operating characteristic (ROC) curves were generated by Gene Expression Profile Interaction Analysis database (GEPIA; http://gepia.cancer-pku.cn/) to assess the diagnostic accuracy of hub genes.The Box-Plot in Expression-DIY in Expression Analysis from GEPIA2 (http://gepia2.cancer-pku.cn/#index)was conducted to verify the expression of hub genes in LC.After that, Kaplan-Meier plotter (http://kmplot.com/analysis/)was applied to reveal the association between the expression levels of the key gene and the survival of patients with LC.The relationship between the expression level of hub genes and the infiltration of immune cell in LC microenvironment was analyzed using Tumor Immune Estimation Resource database (https://cistrome.shinyapps.io/timer/).

Quantitative real time polymerase chain reaction (qRT-PCR)
Total RNA was extracted by Invitrogen™ TRIzol™ Reagent (Invitrogen, Carlsbad, CA, USA), and reverse transcribed to cDNA using PCR Amplifier.qRT-PCR was performed using ABI7500 quantitative PCR (Applied Biosystems, Foster City, CA, USA) with the following procedures: 95 °C for 10 min for initial denaturation, 95 °C for 10 s for denaturation, 60 °C for 30 s for annealing, and 36 cycles in total.GAPDH was used as an internal reference gene.The Ct values were analyzed by the 2 −ΔΔCt method.The experiments were repeated 3 times.Primer (Takara, Otsu, Shiga, Japan) sequence was shown in Table S1.

Western blot
Protein was extracted from cells using lysis buffer (Beyotime), and the protein concentration was determined by bicinchoninic acid quantitative kit (Thermo Fisher Scientific, Waltham, MA, USA).The protein samples were denatured at 95 °C for 5 min.The samples were then separated on the sodium dodecyl sulfate-polyacrylamide gel electrophoresis and were transferred to the polyvinylidene fluoride membranes.After being blocked with 5% skim milk, the proteins were then incubated with anti-BNIPL (1:2000; Abcam, Cambridge, UK) and anti-GAPDH antibody (1:1000, Abcam) at 4 °C overnight.Subsequently, the membranes were incubated with horseradish peroxidase-labeled goat anti-rabbit IgG (1:5000; Abcam) and then visualized using electrochemiluminescence kit (Beyotime).The results were analyzed using an imaging system (Tanon 5200, Shanghai, China).

Proliferation detection
Cell vitality was assessed using cell counting kit-8 (CCK-8) (HY-K0301; MedChemExpress, New Jersey, USA).The cells were plated in 96-well plates (1 × 10 4 /well).CCK-8 assay was performed according to the protocol of instructions.The absorbance at 450 nm was then detected under a microplate reader.

Wound-healing assay
The cells were plated into 6-well plates and scratched using a sterile pipette tip.After washing 3 times with phosphate buffer solution (Thermo Fisher Scientific), cells were added into serum-free medium and incubated in a 5% CO 2 incubator at 37℃.Images of cells were taken at 0 and 24 h using an inverted microscope and analyzed by ImageJ software to calculate the healing rates of cells in both groups, following the formular: healing ratio (%) = (scratch spacing at 0 h -scratch spacing at 24 h)/ scratch spacing at 0 h × 100%.

Transwell assay
Cells (cell density 2 × 10 4 /well) were seeded into the top chamber of transwell insert coated with Matrigel (Corning Inc, New York, USA), and 600 µL 20% fetal bovine serum (Thermo Fisher Scientific) medium was added to the bottom chamber for incubation.After 24 h, cells in the lower surface were fixed in 4% paraformaldehyde (Sigma Aldrich, St. Louis, MO, USA) for 30 min and stained with 0.1% crystal violet (Sigma Aldrich) for 20 min.After washing 3 times with phosphate buffer solution, cells were observed and counted using an inverted microscope.

Statistical analysis
All statistical analysis were conducted using GraphPad Prism 7.0.Data were all represented as mean ± standard deviation.The student's t-test was used to compare the differences between the two groups.P < 0.05 was considered statistically significant.

Identification of DEGs
Total 174 DEGs were screened from GSE143224 dataset, of which 47 genes were upregulated and 127 genes were downregulated.And 575 DEGs were identified from GSE84957 dataset, of which 326 genes were up-regulated and 249 genes were down-regulated.Cluster analysis of the DEGs in the two datasets was carried out to obtain volcano maps, which showed the distribution of DEGs in tumor group compared with the control group (Fig. 1A  and C).The heat map suggested that the samples were clustered with high confidence (Fig. 1B and D).Venn diagram revealed that there were 96 common DEGs in the two datasets (Fig. 2).The data correction results revealed that 25 samples in GSE143224 and 18 samples in GSE84957 were obtained (Figure . S1).The top 25 DEGs (upregulated and downregulated) of GSE143224 and GSE84957 datasets were screened out and displayed in Table S2 and S3.
The GO enrichment results that the hub genes enriched in were mainly "extracellular matrix" and "negative regulation of cell proliferation" (Fig. 5A).PCA analysis showed that the variance interpretation rate of PC1 and PC2 was 87.2%, indicating that hub genes could distinguish between the LC samples from control samples.The scatter plot suggested that the two groups had a good separation, further confirming the validity of the two principal components (Fig. 5B).The expression levels of hub genes in the original sample of GSE84957 dataset were shown in the ridgeline plot (Fig. 5C).

Analysis of BNIPL in LC
Studies have widely reported that KRT4, IGFBP3, MMP10, MMP3, and TGFBI are involved in the pathogenesis of LC and other cancers [15][16][17][18][19].There are few reports focusing on the role of BNIPL in cancers.In present study, we selected BNIPL as the target gene for further analysis.Kaplan-Meier survival curve reflected significant differences of survival in patients with high and low expression (Figure . S2). BNIPL expression was significantly correlated with the infiltration of the immune cells in LC (Figure . S3).Low expression level of BNIPL had a significant association with CD8+, CD4+, natural killer cell, and Macrophage, which promoted the immune infiltration (Fig. 8A-D).

Validation of the hub genes expression
qRT-PCR results suggested that the mRNA expression levels of BNIPL and KRT4 in TU212 and TU686 cells were significantly lower than those in NLECs, while the levels of IGFBP3, MMP10, MMP3 and TGFBI were Fig. 4 Protein-protein interaction network and the identification of hub genes.(A) The protein-protein interaction network of DEGs.(B) The key module where BNIPL resides.The lines between the nodes represent the interactions between genes significantly higher than those in NLECs, which was consist with the bioinformatics analysis results (Fig. 9).

BNIPL inhibits viability, migration, and invasion of LC cells
The results of transfection efficiency suggested the expression levels of BNIPL in TU212 and TU686 cell lines of pcDNA3.1-BNIPL group were significantly increased, compared to the pcDNA3.1-NCgroup (P < 0.001) (Fig. 10A).The CCK-8 results suggested that compared with the pcDNA3.1-NCgroup, the cell viability of TU212 and TU686 cells in pcDNA3.1-BNIPL group was significantly decreased (Fig. 10B).Additionally, the transwell and wound-healing assays revealed that the migration and invasion of TU212 and TU686 cells in the pcDNA3.1-BNIPLgroup were significantly suppressed, compared with the pcDNA3.1-NCgroup (Fig. 10C and  D).

Discussion
Poor prognosis induced by invasion and metastasis of tumor cells is the prevalent issues of advanced and metastatic LC [20].Herein, we identified 6 hub genes (BNIPL, KRT4, IGFBP3, MMP10, MMP3, and TGFBI) of LC with high predictive values using bioinformatics analysis.BNIPL was selected as the key gene, which was correlated with immune cells and promoted the infiltration of immune cell.Moreover, BNIPL was downregulated in LC cells, and overexpression of BINPL inhibited the proliferation, migration, and invasion of LC cells.
Many biomarkers involved in the regulation of malignant tumors have been screened out in previous studies.KRT4 is identified as a hub gene of head and neck squamous cell carcinoma (HNSCC) using CytoHubba after conducting weighted gene co-expression network analysis and differential gene expression analysis between HNSCC and normal tissues [21].IGFBP3 is a regulated gene of p53 tumor suppressor [22].IGFBP3 is identified as a biomarker with significant diagnostic value for early esophagogastric junction adenocarcinoma after ROC curve analysis and consistency index evaluation of the prediction value of line chart [23].After analyzing the data related to oral cancer using the Cancer Genome Atlas, Oncomine, and Kaplan-Meier mapper, the MMP3 and MMP10 are identified as highly altered hub genes of oral cancer [24].TGFBI is identified as one of the hub genes for potential prognostic relevance of breast cancer by conducting weighted gene co-expression network analysis and univariate Cox regression [25].In present research, BNIPL, IGFBP3, MMP10, MMP3, and TGFBI was identified as the hub genes of LC from PPI networks.Overall, BNIPL, IGFBP3, MMP10, MMP3, and TGFBI play crucial roles in the development of LC.
KRT4 inhibits the development of oral squamous cell carcinoma and the expression of KRT4 is downregulated due to m6A methylation at the exon-intron boundary preventing intron splicing of KRT4 pre-mRNA in oral squamous cell carcinoma [15].In patients with LC, serum IGFBP3 levels are significantly higher in the patient group and expresses high diagnostic sensitivity (81.4%) and specificity (80%) [26].MMP1, MMP3, and MMP10 are highly expressed in HNSCC, compared to other cancers by analyzing the data of Oncomine and GEPIA databases.Our findings suggested BNIPL, IGFBP3, MMP10, MMP3, and TGFBI had high sensitivity and feasibility It has been reported that high levels of T cells infiltration in the tumor microenvironment (TME) have diagnostic value for cancers, which plays a crucial role in tumor cells proliferation, invasion, and migration [27,28].Study has reported that infiltration of immune cell in TME is related to clinical prognosis of LC patients: response to induction chemotherapy is better due to less macrophage infiltration and more T cell infiltration in TME [29].In LC, TNNT3, TNNI2, Desmin, matrix metallopeptidase 9, and cytotoxic T lymphocyte antigen 4 are positively correlated with macrophages and dendritic cells, while negatively correlated with CD4 + T cells [30].In our research, the low expression of BNIPL had a significant correlation with the expression of T cells (CD8 + and CD4+) and was correlated with other immune cells in TME.In brief, BNIPL can promote infiltration of immune cell in the TME of LC.
Several studies investigate the effect of BNIPL-related molecules on cancer cell function: growth, migration, and invasion by regulating expression levels.In Hep3B cells, BNIPL-1 upregulates the expression of p16INK4, interleukin-12, TRAIL, and lymphotoxin β genes involved in growth inhibition or apoptosis, and downregulates PTEN expression involved in cell proliferation, which suggests that BNIPL-1 may inhibit Hep3B cells growth through cell cycle arrest or apoptotic cell death pathways [31].In BNIPL-2 transfected Hep3B-Tet-on cells, there are 8 genes involved in cell apoptosis and growth inhibition up-regulating and 7 genes involved in cell proliferation down-regulating, indicating that BNIPL-2 induces apoptosis by regulating the expression of genes involved in apoptosis, growth inhibition, and cell proliferation [32].Poor prognosis induced by ERBB2 can alter the BNIPL expression, causing an anti-apoptotic phenotype [33].In addition, overexpression of CD44 restores the cell proliferation suppressed by BNIPL-2 and BNIPL-2, which can promote migration, invasion, and metastasis of colorectal cancer cells via CD44 [34].In our research, overexpression of BNIPL inhibited the proliferation, migration, and invasion of LC cells.Study has concluded that BNIPL-2 may be a linker protein located at the front of the Bcl-2 pathway for DNA fragmentation and Cdc42 signaling, and for morphological changes during apoptosis [35].This inspires us that BNIPL may play an important role in cancer cells in two pathways: regulating DNA breakage and vesicle formation in apoptotic cells.In addition, the interaction of BNIPL with both MIF and GFER proteins maintains homeostasis between cell proliferation and apoptosis [36].All in all, BNPL exerts an inhibitory effect on the progression of LC.
Taken together, we identified 6 hub genes based on bioinformatics analysis, which had diagnostic and predictive values of LC.In addition, we demonstrated that BNIPL could promote immune cells infiltration in LC and had an inhibitory effect on the proliferation, migration, and invasion of LC cells.Our findings indicate that BNIPL is a promising molecular target and a potential therapeutic for LC, which can provide a new insight for diagnosis and targeted therapy of LC.

Fig. 2 Fig. 1
Fig. 2 Venn diagrams of common DEGs screened from the GSE143224 and GSE84957 datasets.A total of 96 overlapping DEGs are shown in the crossing section

Fig. 3
Fig. 3 Enrichment analysis of common DEGs.(A) Gene Ontology (GO) enrichment bar chart.The abscissa is GO term and the ordinate is -log10 (p-value) of enrichment in each term.(B) Kyoto Encyclopedia of Genes and Genomes (KEGG) bar chart; (C) GO enrichment analysis bubble map, color depth of each node represents corrected p value, the size of the node refers to the number of genes involved; (D) KEGG bubble diagram

Fig. 5
Fig. 5 Analysis of hub genes.(A) GO enrichment chord map, including three parts: logFC is fold-change of genes; other columns are the GO term, and the different links of the gene indicate whether the gene is in this GO term.(B) Principal component analysis of hub genes.PC1 and PC2 on the axes in the figure are the first and second principal components (the variance explanation rate of differences by potential variables); Dots represent samples, and different colors represent different groups.(C) Expression of ridgeline plot of hub gene.The horizontal coordinate is the gene expression, the shape of the mountain represents the dispersion between a set of data, and the height is the number of samples corresponding to the gene expression

Fig. 6
Fig. 6 Receiver operating characteristic curves of hub genes.(A) ROC curves are plotted using expression of BNIPL in GSE84957dataset.(B) ROC curves are plotted using the expression of KRT4 in GSE84957 dataset.(C) ROC curves are plotted using the expression of IGFBP3 in GSE84957 dataset.(D) ROC curves are plotted using the expression of MMP10 in GSE84957 dataset.(E) ROC curves are plotted using the expression of MMP3 in GSE84957 dataset.(F) ROC curves are plotted using the expression of TGFB1 in GSE84957 dataset.The horizontal coordinate is false positive rate and the vertical coordinate is true positive rate, indicating the expression level of the gene in the sample data